地図データの結合

maptools パッケージの spRbind 関数を用いてshapeデータの結合する。 その際、spChFIDs 関数で固有のIDをデータに割り振る必要がある。
例として、統計局の 地図で見る統計(統計GIS) からダウンロードした広島市中区・東区・南区・西区の世界測地系の地図データを結合する。

library(maptools)

# 世界測地系
pj <- CRS("+proj=longlat +datum=WGS84") 

# ダウンロードしたshapeファイルのリスト
filelist <- c(
  "map/h22ka34101.shp", # 広島市中区
  "map/h22ka34102.shp", # 広島市東区
  "map/h22ka34103.shp", # 広島市南区
  "map/h22ka34104.shp"  # 広島市西区
  )

# 1つの地図データに結合
d <- NULL
for(filename in filelist){
  tmp <- readShapeSpatial(filename,proj4string=pj)  # 町丁・字境界データ
  if(is.null(d)) num.d <- 0 else num.d <- nrow(d)
  tmp <- spChFIDs(tmp,as.character(num.d+1:nrow(tmp@data)))
  if(is.null(d)) d <- tmp else d <- spRbind(d,tmp)
}
d <- subset(d, d@data$HCODE == 8101)  # 水上の境界を除外

# 地図を描く
plot(d)

plot of chunk unnamed-chunk-1